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SUMMARY 


A method for generating three-dimensional, finite-difference grids about compli- 
cated geometries by using Poisson equations is developed. The inhomogeneous terms 
are automatically chosen such that orthogonality and spacing restrictions at the body 
surface are satisfied. Spherical variables are used to avoid the axis singularity, 
and an alternating-direction-implicit (ADI) solution scheme is used to accelerate the 
computations. Computed results are presented that show the capability of the present 
method. Since most of the results presented in this paper have been used as grids 
for flow-field computations, this is indicative that the present method is a useful 
tool for generating three-dimensional grids about complicated geometries. 


INTRODUCTION 


In recent years, the development of techniques for generating curvilinear coor- 
dinates has received considerable attention because one of the most important require- 
ments in making finite-difference computations of a flow field involves properly 
locating the nodal points in the flow field (ref. 1). For two-dimensional-flow com- 
putations, highly sophisticated grids have been obtained through a variety of grid- 
generation techniques, such as conformal mappings, algebraic constructions, and 
partial differential equation solutions. Among these, grid-generation schemes based 
on solution of the Poisson equations are popular and have been widely used 
(refs. 2-4). In this method, the coordinates of grid points in the physical domain 
are computed as a solution of a set of Poisson equations. The source terms appearing 
in the equations, which work as forcing functions, control the grid size or shape or 
both near the boundaries. Well-formulated automatic grid-control techniques have 
been developed for two-dimensional problems, and versatile program codes, such as 
NASA’s GRAPE code (ref. 4) developed by Sorenson, are available and have been used 
successfully for many flow-field computations (refs. 5-7). 

The development of grid-generation techniques for three dimensions somewhat lags 
the development of techniques for two dimensions. However, the recent rapid progress 
of high-speed computers makes it possible to perform even three-dimensional Navier- 
Stokes computations, and, hence, three-dimensional grid generation has become the key 
factor for flow-field computations over complex geometries. During the past few 
years, several papers have been written on this subject (refs. 8-10). So far as the 
elliptic partial differential grid-generation techniques are concerned, the essential 
idea is the same as that of two dimensions, although the actual formulation is much 
more complicated. 

The present paper adopts the idea of the two-dimensional code GRAPE and extends 
it to three dimensions. Since the warped spherical coordinates are adequate for the 
three-dimensional flow-field computations, the focus is on constructing such grids. 
There are two unique features in the present paper: the use of the spherical coordi- 

nate variables as independent unknowns and the use of the alternating-direction- 
implicit (ADI) scheme as a solution method. The use of the spherical variables 


removes the axis-singularity problem. The efficiency of the ADI scheme in solving 
the Laplacian-type equations is well known. The constraints of the surface orthogonal- 
ity and the arc-length control are imposed as source terms appearing in the right- 
hand side of the equations to be solved. Details of the formulations and the solution, 
as well as computed grid examples, are presented in the following sections. 

The author wishes to express his gratitude to Mr. J. 0. Bridgeman, who furnished 
his ADI solution code, and to Dr. P. Kutler and Mr. R. L. Sorenson for their valuable 
suggestions during the course of this study. 
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BASIC EQUATIONS AND SOLUTION PROCEDURE 


Although a warped spherical grid is one of the most efficient systems from a 
computational viewpoint, generating such a grid using partial differential equations 
is made difficult because of the trouble associated with the axis singularity. Steger 
(private communication) pointed out that one way to overcome this problem is to solve 
the Cartesian-type Poisson equations with spherical unknown variables. According to 
his suggestion, the mapping between the physical space (p»0,<|>) and the computational 
space (£»U,C) satisfies the following equations (see fig. 1): 

5 pp + «ee + ^ = 
n PP + n ee + n H = 

? pp + ? ee + 4 = 
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The left-hand sides do not correspond to a true Laplacian operator in spherical coor- 
dinates. Equation (1), however, satisfies the maximum principle and generates a 
smoothly varying grid. Interchanging the dependent and independent variables in 
equation (1) results in the following nonlinear equation: 
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Source terms P, Q, and R which appear in the right-hand side of equation (2) 
are determined during the solution process such that the converged solution satisfies 
(1) the orthogonality of the grid at the body surface and (2) the specified grid 
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spacing adjacent to the body in the ^-direction. The conditions for grid orthogonal- 
ity at the body surface are 
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These two conditions indicate that the vector normal to the £ = constant plane 

and the vector normal to n ~ constant plane are orthogonal to the vector normal to 

£ = constant plane. The orthogonality conditions used here are in a different form 

from those of reference 10, but the results are identical. The third condition for 
determining the source terms is control of the grid spacing: 



(3c) 


where S s is the arc length of the ^-direction grid adjacent to the body surface, 
The procedure for determining P, Q, and R from these conditions is described as 
follows. 


Equations (3a) and (3b) are rewittten as 
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By using 0^, <j>£, and p^-, these equations are re-expressed as 
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The relationship between the (p , 0 , <p) coordinates and the (x,y,z) coordinates intro- 
duces the following relation: 

x = (sin 0 cos (f>)p + (p cos 0 cos <J>)0 + (-p sin 0 sin </>)</> 

b b b b 

y^ = (sin 0 sin <j>)p^ + (p cos 0 sin <p)d^ + (p sin 0 cos <j>)<|> ^ (6) 

z = (cos 0)p + (-p sin 0)0 
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Thus, the third constraint, equation (3c), can be expressed as 
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Now, equations (5a), (5b), and (7) can be solved for p^, 0^, and Practically, 

solving equations (5a) and (5b) for 0^ and <j>£ and inserting them into equation (7) 
results in a quadratic equation for p^, which is easy to solve, 1 Since every coeffi- 
cient depends on the body shape only, these equations can be solved a priori and 
desirable p^, 0^, and <J>£ are obtained before starting the solution process of 
equation (2) . 

The source terms are assumed to have the following forms: 

P = P^S.rOe" 31 ? 

Q = Q 1 (5,ri)e' bl? 

R = R 1 (5,n)e _ClC 

Equation (4) indicates that P, Q, and R decay exponentially in the ^-direction. 
Parameters a l9 b l9 and c x control these decay rates. By assuming equation (4), 
source terms in the computational volume are all determined, once source terms at the 
body surface (P 15 Q 19 and R x ) are determined. 
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Solving equation (2) for P, Q, and R results in the following equations; 
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By evaluating equation (5) at the body surface along with equation (4), the following 
relations can be deduced: 


1 Since only one solution of the quadratic equation has meaning, it is necessary 

to choose the appropriate one. For all the cases presented in this paper, the solution 

with a plus sign is chosen. However, it should be carefully chosen in general. 
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Now, the problem to determine P, Q, and R has been reduced to the problem of evalu- 
ating the right-hand side of equation (6) so that three constraints are satisfied. 

Since the £- and q-derivatives at the body surface can be specified from the 
body geometry, only the terms including ^-derivatives a re unknowns. They are p^, 
0 r , <j> r , p ', 0 rr , and <J> rr . Second derivatives are evaluated by the following expres- 
sions (see ref, 4): 
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where subscripts 1, 2, and 3 indicate the indices in the ^-direction (see fig, 2). 
If the first terms in equation (1) are evaluated in the previous solution step, then 
the unknowns to be determined have been reduced to p^, 0^, and <f>£, which can be 
specified in advance from the constraints described above. 


Since instabilities can result if the specified values of p^, 0^, and <j>£ are 
used, the computed values of P, Q, and R are damped in the same way as in ref, 4. 

In each solution step of the ADI scheme, P (for instance) is evaluated as 

p ( n + 1 ) = p (n) + sxgn jMIN[a) p |P - P^ | , P Lim • MAX( | P n | , 1) ] , P - P^j (12) 

where the superscript indicates the number of the iteration and where a)p is the 
underrelaxation parameter. The same procedure is used for Q and R. 


Equation (2) is solved by the ADI scheme with the source-term evaluation process 
shown above. Since the ADI scheme is widely used as a solution method for elliptic 
equations, the details are not shown here. A false time-step is introduced and the 
solution of the original elliptic equation is obtained as the time-asymptotic steady- 
state solution. 
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APPLICATION TO THE GENERATION OF WARPED SPHERICAL COORDINATES 


The essential idea of the transformation between the physical and computational 
space is shown schematically in figure 1. Bilateral symmetry is imposed, so that the 
grid is generated in a half-volume. At the axis, on the body boundary, and on the 
outer boundary, p, 0, and <f> values are specified. At the outflow boundary, p, 0, 
and <J> are extrapolated so that 3x/9n - 0 and 9y/9ry = 0 are satisfied in specified 
z = constant plane. 

, Since convergence difficulties arose in the ADI solution process, the following 
optional modifications were included in the program code. First, time-steps can be 
varied for each of the three equations in equation ( 3 ), for small time-steps are 
required for the p-equation in the initial stage of the solution process. Second, 
the p-transf ormation can be introduced to overcome the trouble appearing near the 
origin when the origin is chosen to be near the body surface. The orthogonality con- 
dition is distorted to some extent when this transformation is introduced. Third, 
the grid can be redistributed in the ^-direction. Because stable computations for 
the very small spacing near the body were difficult, the computations are done for the 
relatively larger spacing near the body surface, and only grid lines in the 
^-direction are kept. Then, the grid points are redistributed from the body surface 
to the outer boundary along the obtained 5 -lines. This process makes it possible to 
obtain Euler and thin-layer Navier-Stokes grids by a single run of the grid- 
generation program. 

It should be noticed that the location of the origin changes the convergence 
history very much in the present method, using spherical variables, whereas the origin 
location does not matter when the Cartesian variables are used. 

With all these modifications, small time-steps had to be used to avoid the insta- 
bilities, and, as a result, convergence was slow and required premature termination of 
the iteration process before the convergence criterion was satisfied. Thus, the grids 
obtained are not strictly normal to the body surface but are adequate for the flow- 
field computations, as is seen in the illustrative examples. 


ILLUSTRATIVE EXAMPLES 


application of this program code to a wide variety of body shapes 
The grid over the body and the outer boundary are given as 
and the initial guess for the inner region grid is straight lines 
two boundaries in all the cases presented here. 


Hemisphere-Cylinder 

The first example is the application of the code to the hemisphere-cylinder. 
Since the body is axisymmetric , the grid generation is rather easy in this example. 

In figure 3, the body shape, as well as the computed grid distribution, is presented. 
The number of grid points is 24 along the axis, 21 in the circumferential direction, 
and 21 in the radial direction. Stretching of the body-grid spacing is introduced in 
the rrdirection so that the grid obtained can be practically used. Since <j) values 
are the same for each E, = constant plane, it is not necessary to solve the 


Examples of the 
are presented below, 
boundary conditions, 
for connecting these 
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<j>-equation, even though the equation is actually solved together. This figure 
includes J = 2, J = 20, and K = 24 grid surfaces, as well as the body grid. 

The redistributed grid, adequate for viscous computations, is shown in fig- 
ure 4(a), and the blown-up view of £ -plane grid is shown in figure 4(b). The number 
of the grid surface in the radial direction is changed from 21 to 31 in this viscous 
grid. The result indicates that the grid obtained is smooth, and the orthogonality 
is almost satisfied in this example. 


Delta Wing without Trailing Edge 

The body shape (fig. 5) is composed of the wing region and the following slab 
region, which is attached to avoid the expected difficulty associated with the wake 
treatment. The number of grid points in this example is 31 in the chordwise direc- 
tion, 31 in the circumferential direction, and 21 from the body to the outer boundary. 

Figure 6 shows the overall view of the grid distribution where the J = 2, 

J = 16, and K = 31 planes are shown, as well as the body grid. Since the focus is 
on the wing region, a very large stretching ratio for the slab region is used in the 
chordwise direction. The grid distribution is obtained only for the forward part of 
the body, and the following grid is obtained by extrapolation, keeping only y and z 

data from the last computed grid plane in the n-direction. This is done in order to 

save computer time and avoid instability. 

The redistributed viscous grid that was actually used for the thin-layer Navier- 
Stokes computations (ref. 11) for the leading-edge-separation vortex study is shown 
in figure 7(a). Close-up views from the front and the top are also shown in fig- 
ures 7(b) and 7(c). Figure 7(b) is the top view of J = 2, J = 16 grid surfaces and 

body grid. Figure 7(c), the view from the front, shows the J = 2, J = 30 , and 
K = 15 surfaces, as well as the body grid. It should be remembered that the K = 15 
crossflow plane is not a true plane but a curved surface bending toward the front, as 
is seen in figure 7(a). The origin was taken to be the apex of the body, and thus 
the weak p-trans formation was introduced in this example. These figures indicate 
that the grid is not strictly orthogonal, but that it is adequate for the flow-field 
computations. The grid distribution caused no troubles throughout the flow-field com- 
putations. Much clustering and stretching is used on the boundary in both the £- and 
ri-directions to make the grid useful for the flow-field computation. Thus, the com- 
putation of the grid distribution in this case is not as easy as the computation for 
the constantly distributed body-grid case. 


Space Shuttle Geometry 

The present method was applied to the complicated geometry of the Space Shuttle. 
The body shape is given by the geometry package called "GE0M3." Complete details of 
the geometry package used here are given in reference 12. The body grid is obtained 
by cubic-spline interpolation, so that the grid is clustered near the wing-tip region, 
etc. The number of grid points is 38 in the chordwise direction, 51 in the circum- 
ferential direction, and 31 from the body to the outer boundary. This is almost the 
maximum number of grid points used for the flow-field computations, using the 
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vectorized version of the thin-layer Navier-Stokes solution code. 2 Because of this 
limitation, the wing region does not have enough points in the chordwise direction. 

A grid having 46 points in the chordwise direction was also generated and showed a 
very similar result. 

Figure 8 shows the geometry and the body grid. The distribution of the grid 
points over the geometry is performed by using an interactive grid-distribution pro- 
gram including a display capability. Figure 9(a) shows the overall view of the 
viscous grid, where j = 2, j = 50, and k = 38 planes are seen. Figures 10(a)-10(d) 
show the view of the crossflow grid from the front. Again, it should be remembered 
that these are not planes but curved surfaces. In each figure, the grid orthogonality 
is not satisfied but the grid is acceptable. The grid distribution is smooth, even in 
the region of wing-body junction where grid generation might be very difficult. 


CONCLUSIONS 


A grid-generation program code has been developed for three-dimensional geome- 
tries. In this code, spherical variables were used to avoid the axis-singularity 
problem, and a set of three-dimensional Poisson equations was solved by the ADI 
scheme. The orthogonality condition and the minimum spacing restriction on the body 
surface were imposed as source terms in the equations. The program code was written 
for generating the warped spherical coordinate system, but can be easily extended to 
other topologies. 

The program code was applied to several geometries, and the results indicate that 
it can generate a grid that is adequate for the flow-field computations, even though 
the strict orthogonality was not satisfied. 


2 The original program code was capable of using only half this number of points. 
A modification of the program was made by the present author, and this number of grid 
points for the flow-field computations became possible using the Solid State Disk 
Device (SSD) on the CRAY 1-S. 
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APPENDIX 


APPLICATION OF THE PRESENT METHOD TO THE THIN, DELTA WING 


The present code was applied to generate the warped spherical coordinate system 
over a thin, delta wing. Since generating thin, sometimes zero-thickness, body geome- 
try is very difficult, the following procedure is used. First, the grid distribution 
is done for a thicker wing having no trailing-edge. The grid-generation process is 
the same as that used in the delta wing example above. The body grid for the thin 
wing, which has zero thickness in the wake region, is distributed so that the spacing 
ratio in both the £- and p-directions are the same as for the thicker wing. Then, 
each 5-line for the thicker body is extended to the corresponding thin-body grid; 
5-lines are smoothly redistributed, using cubic-spline interpolation. A typical 
crossflow plane grid in the wing region is shown in figure 11; figure 12 shows the 
grid in the wake region. Since the body is very thin, the orthogonality is not much 
distorted by the process above. The flow-field computations were actually done and 
were successful. In this topology, the transformation Jacobian becomes singular along 
the line starting at the wing tip in the wake region. This was overcome by special 
treatment of that line. 
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Figure 1.- Correspondence between physical space and computational space. 
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(a) Redistributed viscous grid over hemisphere-cylinder. 
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(b) Close-up view of £ -plane grid. 

Figure 4.- Grid over hemisphere-cylinder for viscous flow computations. 
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Figure 5.- Body geometry. 
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Figure 6.- Overall view of the grid. 
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(a) Overall view of the redistributed viscous grid. 


(b) Close-up view of the grid from the top. 

Figure 7.- Grid over delta wing for viscous flow computations. 
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(c) Close-up view of the grid from the front. 
Figure 7.- Concluded. 
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Figure 9 



Figure 8.- Shuttle geometry and the body grid. 



, - Overall view of the grid (redistributed for the viscous computations). 
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(b) K = 30. 


Figure 10.- Close-up view of the crossflow plane grid. 
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(d) K = 38. 


Figure 10.- Concluded. 
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Figure 11.- Crossflow surface grid for thin wing (wing region). 



Figure 12.- Crossflow surface grid for thin wing (wake region). 
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